#include "gtest/gtest.h"

#include "sgCPardiso.h"

// #include "DataStructure/Common/include/Constant.h"

//#include <complex.h>
#include <complex>

using namespace SG::Algebra;  // Open the namespace
using namespace std::complex_literals;

class CPardisoTest : public ::testing::Test
{
  public:
    static void SetUpTestSuite ()
    {}
};

/*
 * Test the A*x = b given by SuperLU User Guide.
 * where A is a 3X3 matrix and b is a 3X1 matrix.
 *
 * A =
 *   [
 *       06.00,  07.00,  05.00;
 *       07.00,  13.00,  08.00;
 *       05.00,  08.00,  06.00
 *   ]
 *
 * b = [9.0; 10.0; 9.0]
 *
 * OctaveOnline give the solution as followings,
 * x = inv(A)*b = [1; -1; 2]
 *
 * see https://octave-online.net/
 */
TEST_F (CPardisoTest, realTest)
{
    const int N = 3;

    Real_t coeffs[] = { 6.0, 7.0, 5.0, 7.0, 13.0, 8.0, 5.0, 8.0, 6.0 };
    Real_t rhs[]    = { 9.0, 10.0, 9.0 };

    SpMat_t   A (N, N);
    const int row_start = A.row_start ();
    const int row_end   = A.row_end ();

    for (int i = row_start; i < row_end; ++i)
    {
        for (int j = 0; j < 3; ++j)
        {
            A.set (i, j, coeffs[3 * i + j]);
        }
    }
    A.assemble ();

    CPardiso<Real_t> solver;
    solver.compute (A);

    // Check
    auto info1 = solver.info ();
    ASSERT_TRUE (Info_t::Success == info1);

    Vec_t b (N);

    for (int i = row_start; i < row_end; ++i)
    {
        b.set (i, rhs[i]);
    }
    b.assemble ();

    Vec_t x (N);
    solver.solve (b, x);

    // Check
    auto info2 = solver.info ();
    ASSERT_TRUE (Info_t::Success == info2);

    // solution from OctaveOnline
    // x = inv(A)*b
    // see https://octave-online.net/
    Real_t exact[] = { 1, -1, 2 };

    for (int i = row_start; i < row_end; ++i)
    {
        ASSERT_NEAR (exact[i], x (i), TOLERANCE);
    }
}

TEST_F (CPardisoTest, complexTest)
{
    const int N = 3;

    Complex_t coeffs[] = { 6.0I, 7.0I, 5.0I, 7.0I, 13.0I, 8.0I, 5.0I, 8.0I, 6.0I };
    Complex_t rhs[]    = { 9.0I, 10.0I, 9.0I };

    ComplexSpMat_t A (N, N);
    const int      row_start = A.row_start ();
    const int      row_end   = A.row_end ();

    for (int i = row_start; i < row_end; ++i)
    {
        for (int j = 0; j < 3; ++j)
        {
            A.set (i, j, coeffs[3 * i + j]);
        }
    }
    A.assemble ();

    CPardiso<Complex_t> solver;
    solver.compute (A);

    // Check
    auto info1 = solver.info ();
    ASSERT_TRUE (Info_t::Success == info1);

    ComplexVec_t b (N);

    for (int i = row_start; i < row_end; ++i)
    {
        b.set (i, rhs[i]);
    }
    b.assemble ();

    ComplexVec_t x (N);
    solver.solve (b, x);

    // Check
    auto info2 = solver.info ();
    ASSERT_TRUE (Info_t::Success == info2);

    // solution from OctaveOnline
    // x = inv(A)*b
    // see https://octave-online.net/
    Complex_t exact[] = { 1, -1, 2 };

    for (int i = row_start; i < row_end; ++i)
    {
        ASSERT_NEAR (exact[i].real (), x (i).real (), TOLERANCE);
        ASSERT_NEAR (exact[i].imag (), x (i).imag (), TOLERANCE);
    }
}